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HIGHLIGHTS 


►  We  extended  the  physics-based  single  particle  model  for  higher  charge-discharge  rates  up  to  5C. 

►  The  improved  single  particle  (ISP)  model  was  compared  with  full-order  pseudo  two  dimensional  (P2D)  model. 

►  The  ISP  model  could  predict  the  P2D  model  results  up  to  5C  while  the  simulation  time  was  reduced  by  a  factor  about  5. 
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An  isothermal  reduced  order  physics-based  Li  ion  battery  model  is  developed  by  incorporating  solution 
phase  charge  and  material  balances  into  single  particle  (SP)  model.  Li  ion  concentration  and  potential 
profiles  in  electrolyte  phase  are  approximated  by  polynomial  functions.  A  cubic  polynomial  is  used  for 
electrolyte  concentration  and  potential  inside  electrodes,  while  separator  liquid  phase  potential  and 
concentration  are  estimated  by  parabolas.  Diffusion  inside  solid  particles  is  also  simplified  using  an 
approximate  solution  based  on  analytical  solution  for  the  solid  concentration.  Applying  the  aforemen¬ 
tioned  reduction  techniques  decreases  the  degree  of  freedom  of  full  order  model  by  more  than  100,  while 
the  cell  voltage  relative  error  of  proposed  model  is  less  than  1  %  at  different  charge-discharge  rates  up  to  5C. 
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1.  Introduction 

Due  to  high  energy  density  of  Li  ion  cells  with  respect  to  other 
kind  of  batteries  (three  to  four  times  of  Lead  acid  and  Nickel 
cadmium  and  twice  of  Nickel  metal  hydride)  rechargeable  Li  ion 
batteries  are  used  extensively  in  different  markets  and  applications 
ranging  from  consumer  electronics  (e.g.  laptops,  cell  phones)  to 
automotive  (e.g.  Hybrid  Electrical  Vehicles)  and  aerospace  (e.g. 
satellite  power  sources)  [1  ].  However,  Li  ion  cells  are  unforgiving  if 
operated  outside  of  tight  safe  operating  area  [2].  In  order  to  guar¬ 
antee  safe  charge-discharge  operations  of  a  Li  ion  cell,  most  battery 
systems  consist  of  a  battery  management  system  (BMS)  to  control 
charging-discharging  [2].  The  core  of  a  BMS  is  a  battery  model  that 
finds  the  relationship  between  currents  and  voltages  measured  at 
terminals  [3]. 

Li  ion  cell  models  are  categorized  into  two  main  groups.  The  first 
type  is  empirical  approaches  that  rely  on  capacitor-resistor 
networks  such  as  equivalent  circuit  model  (ECM)  [4].  Due  to 
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simplicity  and  low  computational  burden,  empirical  models  have 
been  used  widely  in  online  implementations  [4].  However,  physical 
insight  of  the  cell  is  neglected,  resulting  in  prediction  errors  when 
used  over  a  broad  operating  region  [3].  The  second  type  includes 
physics-based  models  like  Pseudo  Two  Dimensional  (P2D)  [5,6] 
where  distribution  of  essential  states  inside  the  cell  (e.g.  solid  and 
liquid  phase  Li  ion  concentrations  and  potentials)  along  cell  sand¬ 
wich  is  obtained  by  solving  material  and  charge  balances.  Thus, 
physics-based  models  can  be  applied  in  different  operation  modes 
and  also  capacity  fade  phenomena  are  included  in  the  model  without 
any  difficulty  [7].  Nevertheless,  the  complexity  and  computational 
cost  restricts  the  use  of  these  models  for  onboard  applications. 

In  order  to  make  rigorous  physics-based  models  more  compu¬ 
tationally  efficient,  various  reduced  models  have  been  proposed 
[8,9,11-16].  The  most  simplified  model  developed  by  Atlung  et  al. 
[8]  and  later  extended  by  Haran  et  al.  [9]  is  single  particle  (SP) 
model  where  each  electrode  is  represented  by  a  single  spherical 
particle  and  Li  ion  concentration  and  potential  inside  electrolyte 
phase  are  neglected.  Although  the  SP  model  computation  time  is 
comparable  with  empirical  ECM  [10],  it  is  limited  for  low  rate 
operations  (e.g.  <1.0C)  [11].  V.R.  Subramanian  et  al.  approximated 
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microscale  diffusion  equation  in  a  spherical  electrode  particle 
based  on  polynomial  profile  [12].  They  coupled  the  approximated 
solution  with  governing  equations  for  the  macroscale  to  predict  the 
cell  voltage  as  well  as  electrolyte  concentration.  Although  it  was 
shown  that  the  number  of  differential  algebraic  equations  (DAEs) 
reduced  by  order  of  10,  the  resulting  DAE  is  still  large.  In  the  next 
work  by  V.R.  Subramanian  et  al.,  the  full  order  DAE  system  was 
reduced  further  by  approximating  electrolyte  concentrations  inside 
electrodes  and  separator  using  polynomial  functions  [13].  They 
showed  that  approximated  solutions  for  different  quantities  such 
as  electrolyte  concentrations  and  potential  agreed  well  with  full 
order  profiles  under  1C  discharge  rate.  However,  uniform  electro¬ 
lyte  current  densities  were  assumed  in  each  electrode.  As  we  will 
show  later  in  this  paper,  this  assumption  is  true  only  for  low 
charge-discharge  rates. 

A  physics-based  reduced-order  model  was  proposed  by  L.  Cai 
and  R.E.  White  based  on  proper  orthogonal  decomposition  [14]. 
There  were  50  equations  in  the  reduced  model  (compared  to 
14,580  equations  in  high-order  dimensional  discrete  model)  that 
saved  the  computational  time  of  the  rigorous  model  by  a  factor 
about  7.  Smith  et  al.  developed  a  low  order  electrochemical  model 
under  assumptions  of  quasi-linear  behavior  (i.e.,  constant  elec¬ 
trolyte  properties  and  linearized  Butler  Volmer  (BV)  kinetics)  and 
reaction  current  decoupled  from  electrolyte  concentration  [15]. 
They  reduced  a  full  order  impedance  transfer  matrix  to  a  12th 
order  state  variable  model  that  predicted  the  cell  voltage  to  within 
1%  for  pulse  and  constant  current  profiles  at  rates  up  to  50  C. 
However,  the  validity  of  the  reduced  model  was  lost  when  the 
electrolyte  became  depleted  (Ce  <  0.15  Ce,o).  Smith  et  al.’s  model 
was  further  simplified  by  D.  Di  Domenico  et  al.  assuming 
a  constant  electrolyte  concentration  and  neglecting  the  solid 
concentration  distribution  along  each  electrode  (resulting  in 
constant  current  density  in  BV  kinetics)  [16].  Their  reduced  model 
is  similar  but  higher  order  than  SP  model.  Recently,  T.-S.  Dao  et  al. 
presented  a  simplified  model  using  a  volume  average  technique 
for  solid  phase  concentration  and  Galerkin’s  approximation 
for  electrolyte  phase  concentration  and  potential  [17].  The 
current  density  was  also  assumed  to  be  constant  along  each 
electrode.  They  obtained  a  set  of  2 N  +  10  (N  is  number  of  node 
points  for  Galerkin’s  approximation)  DAEs  where  N  =  4  resulted 
a  reduced  model  in  good  agreement  with  rigorous  model  for  1  C 
discharge  rate. 

In  this  paper,  the  isothermal  SP  model  without  capacity  fade 
[10]  is  improved  to  predict  the  cell  voltage  at  higher  charge- 
discharge  rates  up  to  5C  where  the  Li  ion  concentration  in  elec¬ 
trolyte  phase  (near  to  cathode/anode  current  collector  during 
discharge/charge)  depletes  near  to  zero  (Ce  ~  0.001  Ce,o).  Distri¬ 
bution  of  electrolyte  concentration  and  potential  along  electrodes 
and  separator  are  incorporated  with  the  SP  model  equations  by 
using  polynomial  approximations  for  the  electrolyte  concentra¬ 
tions  and  potentials.  The  new  efficient  reduced  model  consists  of 
only  13  DAEs  where  cubic  and  quadratic  polynomials  are  applied 
for  solution  phase  transport  and  conduction  in  electrodes  and 
separator  respectively.  As  mentioned  in  Ref.  [13]  the  proposed 
method  is  slightly  different  from  standard  collocation  procedures 
because  governing  equations  are  also  volume  averaged  to  find 
polynomial  coefficients.  In  order  to  enhance  the  accuracy  of  solid 
phase  surface  concentration,  the  volume  average  technique  used 
in  Ref.  [10]  is  replaced  by  an  approximated  solution  developed  by 
M.  Guo  and  R.E.  White  [18].  Comparison  of  constant  current 
charge-discharge  data  obtained  by  full  order  P2D  with  improved 
single  particle  (ISP)  model  shows  that  there  is  a  good  agreement 
between  the  reduced  and  the  rigorous  models  while  the  compu¬ 
tation  time  of  the  new  model  is  greatly  less  than  the  full  model  (a 
factor  of  5). 


2.  Physics-based  model 

In  this  section,  the  assumptions  and  governing  equations  of  full 
order  P2D  model  for  Li  ion  concentrations  and  potentials  inside 
solid  and  solution  phases  are  first  presented.  Then,  the  underlying 
assumptions  and  restrictions  of  simplified  SP  model  are  addressed. 
Finally,  the  SP  model  is  extended  to  higher  rates  by  coupling  the 
electrolyte  charge  and  material  balances  to  solid  phase  diffusion.  In 
all  models,  the  BV  expression  is  used  to  predict  the  rates  of  Li  ions 
intercalation/deintercalation  reactions. 


2.1.  Pseudo  two  dimensional  model 


Fig.  1  presents  a  schematic  of  a  Li  ion  cell  composing  of  two 
porous  electrodes  (solid  matrix  inside  an  electrolyte  solution)  and 
a  separator  (electrolyte  solution).  Most  physics-based  full  order  Li- 
ion  cell  models  rely  on  the  rigorous  model  of  Doyle  et  al.  [5]  and  T.F. 
Fuller  et  al.  [6].  Transport  of  Li  ions  in  the  electrolyte  phase  is 
considered  only  in  x  direction  (i.e.  from  the  cathode  through  the 
separator  into  the  anode  during  charge  and  vice  versa  during 
discharge)  since  the  lengths  of  the  cell  in  y  and  z  directions  are 
much  greater  than  the  cell  thickness  (>1000).  In  the  solid  phase  Li 
ions  are  assumed  to  react  at  the  surface  and  move  only  in  r  direc¬ 
tion.  Thus,  simulation  of  the  cell  can  be  regarded  as  a  packed  bed 
reactor  model.  In  addition  to  material  balances  of  Li  ions  inside 
liquid  and  solid  phases,  charge  balances  for  each  phase  are  required 
in  order  to  obtain  the  potential  distribution  across  the  cell.  Hence, 
there  are  four  quantities  describing  the  cell  electrochemical 
behavior  including  solid  and  electrolyte  Li  ion  concentration  and 
potential.  These  variables  are  coupled  by  a  charge-transfer  kinetic 
resistance  at  the  particle  surface.  In  the  following,  governing 
equations  and  assumptions  are  presented. 


2.1.1.  Solid  concentration 

Li  ions  travel  inside  spherical  solid  particles  by  diffusion  process 
given  by  the  Fields  second  law: 


9Cs,i 

dt 


Ds,i  d  /20QA 

r2  dr  \  dr  J 


i  =  p,n 


(la) 


with  the  initial  and  boundary  conditions: 


CsJ(r,t  =  0)  =  CS  i  o  0  <  r  <  Rj 
9Cs>i 


Or 

—Dsj 


=  0 


t>  0 


r= 0 

dCsf 

dr 


r=Ri 


=  Ji(X,t) 


(lb) 


where  Ji  is  the  pore  wall  flux  at  the  interface  between  electrolyte 
and  particles.  DS) i,  the  diffusion  coefficient  of  Li  ions  inside  solid 
particles,  is  a  constant. 


2.1.2.  Solid  potential 

The  potential  distribution  in  the  solid  phase  of  each  electrode  is 
obtained  by  Ohm’s  law  together  with  the  boundary  conditions  at 
two  ends  of  the  cell  and  the  electrode/separator  interfaces: 


d  2(k_ 

=  a>FJ>(x,t)  i  =  P,n 

Cathode  :  —  ffeff,p 


=  I. 


Anode :  ^ 

dx 


x=0 

=  0 


app 


x — Lp  +LS 


9</>p 

0X  x—Lp 
Pn\x=L  =  0 


=  0 


(2) 
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where  <reff,i  is  constant  effective  solid  conductivity  and  ax,  specific 
interfacial  area,  is  determined  by: 


9Ce,s  _  1  9  f  n  9Ce,s^\ 

“a F  _  V  e,s_a3T ) 


(3b) 


3(l-  ei  -  £f,i) 


•  Anode: 


where  fc-f.i  is  volume  fraction  of  fillers.  The  charge  flux  at  the  current 
collector/cathode  interface  is  equal  to  the  applied  current  to  the  cell 
while  there  is  no  charge  flux  at  the  electrode/separator  interfaces. 
The  potential  at  the  current  collector/anode  interface  is  arbitrary 
set  to  zero. 

2.1.3.  Electrolyte  concentration 

The  material  balance  on  the  salt  in  the  electrolyte  phase  for 
different  regions  is  obtained  by  assuming  concentrated  solution 
theory  (a  binary  salt  plus  solvent)  [19]: 

•  Cathode: 


9Ce,n  _  1  9  fn  9Ce,n^  (1  f+)<Wn(*?0  ^ 

~dF  ~  £^dx  {De  n~dTJ  +  ^  (3C) 

where  the  transference  number  of  Li  ions,  tp,  is  assumed  to  be 
constant  and  the  effective  diffusion  coefficient  of  the  lithium  salt, 
De> i,  depends  on  the  electrolyte  concentration,  temperature, 
porosity  and  Bruggeman  number  [20] 


De.i  =  f|brugg'D|,  i  =  p,s,n 


Di  =  10 


-4.43- 


- 54 - 0.22xl03CeI(x,t) 

T-229-5.0xl03Cei(x,f)  e,IV  ’  ’ 


(4) 


The  following  boundary  conditions  are  applied  for  the  electro¬ 
lyte  concentration: 


8Ce,p  _10fn  aCe.p^  ,  (l-t°)ap/p(x,f) 
9 1  £pdx\  e,p  9x  )  +  Cp 


•  No  mass  flux  at  two  ends  of  the  cell  in  x-direction 

(3a) 


9Ce,p(x,  t) 
Dep  9x 


•  Separator: 


x=0 


(5a) 
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—Den 


aCe,n(X,  0 


Ox 


=  0 


x=L 


(5b) 


Continuity  of  the  concentration  and  the  flux  at  the  cathode/ 
separator  interface: 


Ce,p(Lp,t)  -  Ce,s(£p,t) 


-D, 


0Ce,p(x,  t) 


e,P“ 


OX 


=  —De 


ace,s(x,  t) 


x — Lp 


Ox 


x — Lp 


The  following  boundary  conditions  are  used  for  the  electrolyte 
potential  inside  the  cell: 

•  No  charge  flux  at  two  end  of  the  cell: 


(5c) 

/  a0e,p 

ke  p  9x 

=  0 

x=0 

(8a) 

(5d) 

i  ^0e,n 

ke  n  9x 

i  o 

x=L 

(8b) 

•  Continuity  of  the  concentration  and  the  flux  at  the  separator/ 
anode  interface: 


•  Continuity  of  potential  and  charge  flux  at  the  cathode/sepa¬ 
rator  interface: 


Ce,s  (Op  +  f) 

—  Ce,n(£p  +  f) 

(5e) 

^e,p(0pO)  —  ^e,s(0pO) 

(Be) 

n  9Ce,s(x,  t) 

n  9Q,n(X,  t) 

(5f) 

i  ^^e,p  i  ^^e,p 

<e  p  0x  v  ,  “  ^  0X  v  , 

e’s  ax 

e,n  ax 

x=Lp+Ls  ax 

x — Lp  +LS 

(8d) 

2 A. 4.  Electrolyte  potential 

The  variation  in  electrolyte  potential  inside  the  cell  is  deter¬ 
mined  by  applying  Ohm’s  law  having  assumed  the  electrolyte  is  an 
ideal  solution: 


•  Continuity  of  the  potential  and  charge  flux  at  the  interface 
separator/anode: 


•  Cathode: 


_a_ 

ax 


+  8  -  f/<e  p  8Ce’p) 

^9xVCe,p  0X  ) 


apFJp(x >  t) 


^e,s  (Lp  +  Ls,t)  —  4>e.n  (Lp  +  Ls,  t) 


(6a) 


-ke 


d<t>e 


dX 


x — Lp+Ls 


—  —ke,n 


^0e,n 


ax 


x — Lp+Ls 


(Be) 

(Bf) 


•  Separator: 


a 

ax 


+  lgA  (hi 9CeA 

M0xVCe,s  9x  J 


=  0 


(6b) 


•  Anode: 


_a_ 

ax 


n  d  //<e,n  aCe,n\ 

PaxVCe, n  ax  J 


an^/n(X,  t) 


where  (3  is  equal  to  2FgT(l  -  t°)/F  and  the  specific  conductivity  is 
a  function  of  the  electrolyte  concentration,  temperature,  porosity 
and  Bruggman  number  [20]: 


2 A. 5.  Pore  wall  flux 

The  BV  kinetic  expression  is  used  to  predict  the  rates  of  the  Li  ion 
deintercalation/intercalation  reactions  for  each  electrode: 

Ji(xi  0  —  2/<;  ^Cs  max  / 

-  \0fY  r  \°'5r0.5  ■  uf°-5F  ^ 

Xs, surf, iQ, max, ij  (^Xs,  surf, iQ, max,  i J  j Vi) 

(9) 

The  electrode  overpotential  for  the  lithium  ion  intercalation/ 
deintercalation  reactions  is  given  by: 

Vi  —  0i  —  0e,i  —  (xs,surf,i)  09) 

where  U[  denotes  electrode’s  open  circuit  potential  [7].  The  cell 
voltage  is  calculated  by  the  following  equation: 


kei  =  £°mgg’I<h  i  =  p,s,  n 

_4  / -10.5  +  0.668  x  10_3Cez  +  0.494  x  lO”6^.  +  0.074F  -  1.78x  \ 

I<j  =  10  Ce  i  ^  10-5q  T  _  8  86  x  i0-ioc2 ’-  .j  _  6  96  x  \q-5T2  +  2  80  x  10-8Ceir2  j 


2 


(7) 
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^cell  —  ^p|x_q  ^nlx=L  (H) 

2 A. 6.  DAE  system 

In  order  to  solve  the  system  of  coupled  PDEs  of  the  full  order 
model,  discretization  along  x  axis  in  the  electrodes  and  the  sepa¬ 
rator  in  association  with  discretization  along  r  axis  for  the  solid 
diffusion  are  required.  The  P2D  model  is  solved  in  COMSOL  Inc. 
Multiphysics  3.5a  by  using  70,  25,  74  elements  along  x  axis  for 
cathode,  separator  and  anode  respectively.  The  particles  in  solid 
phase  are  divided  into  200  elements.  If  the  element  type  is  selected 
as  Lagrange- Quadratic  (default),  the  model  consists  of  2690 
degrees  of  freedom.  The  model  parameters  used  in  this  study  are 
listed  in  Table  1. 

2.2.  Single  particle  model 


dx, 


s,avg  ,i 


dt 


~3 l 

^2  Q,  max,  i 


-JfRi 

*s,surf,i  ^s,avg,i"^  cn  r 

JL/s,i'“-s,max,i 


,1  =  p,n 


(12) 


where  Jj  denotes  the  total  flux  of  Li  ions  intercalate  into  or  dein- 
tercalate  from  an  electrode: 


fT  =  {£PP 
Jp  FSP 

tT  =  ~^pp 

M  FSn 


(13) 


where  Si  is  the  electrode  total  active  surface  area  that  can  be 
determined  as  the  product  of  the  specific  interfacial  area  and  the 
electrode  volume: 


SP  model  considers  each  electrode  as  a  single  spherical  particle 
whose  area  represents  the  porous  electrode  active  surface  area. 
That  is,  it  is  assumed  that  all  of  the  particles  in  an  electrode  behave 
the  same  way  and  that  the  current  being  passed  through  the 
electrode  is  distributed  uniformly  over  all  of  the  particles.  More¬ 
over,  the  variation  of  electrolyte  concentration  and  potential  with 
position  and  time  are  ignored.  Therefore,  the  system  of  full  order 
PDEs  is  reduced  to  two  solid  diffusion  PDEs  (Eq.  (1))  in  which  each 
electrode  regarded  as  spherical  particles.  The  solid  diffusion 
equation  can  also  be  simplified  by  applying  a  volume  average 
technique  where  a  polynomial  is  used  to  approximate  the 
concentration  profile  inside  the  particle  [21].  The  PDE  and 
boundary  conditions  in  Eq.  (1)  are  converted  to  the  following  DAE 
system  using  a  parabolic  polynomial  approximation: 


Table  1 

Pseudo  two  dimensional  model  parameters. 


Parameter 

Value 

Unit 

Ce,o 

1000 

mol  m  3 

kp 

6.667e-ll 

m2-5mora5s_1 

kn 

1.764e-ll 

m2-5  mol-0  5  s— 1 

RP 

8.5e-6 

m 

Rn 

12.5e-6 

m 

Ds,p 

le-11 

m2  s-1 

f^s.n 

5.5e-14 

m2  s-1 

Q,max,p 
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b 
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In 
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m 

t? 
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bruggs 
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bruggn 
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m3 

Vn 
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m3 

Sp 
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m2 

Sn 

0.75 

m2 

T 

298.15 

K 

Kg 

8.3143 

J  mol1  K  1 

F 

96,487 

C  mol-1 

Cell  capacity 

1.656 

Ah 

Cathode  initial  state  of  charge  for  charge 

0.95 

- 

Anode  initial  state  of  charge  for  charge 

0.05 

- 

Cathode  initial  state  of  charge  for  discharge 

0.465 

- 

Anode  initial  state  of  charge  for  discharge 

0.756 

- 

Si  =  OiVi  (14) 

In  summary,  the  solid  surface  concentration  and  potential  are 
functions  of  time  only  and  the  electrolyte  concentration  and 
potential  are  constant  restricting  the  model  usage  for  low  rate 
operation  (<1.0C)  [11].  The  DAEs  of  SP  model  can  also  be  solved 
analytically  when  the  current  is  constant,  the  model  does  not 
include  capacity  fade,  and  cathodic  and  anodic  transfer  coefficients 
are  the  same  (a  =  0.5)  [10]. 

2.3.  Improved  single  particle  model  (ISPM) 

This  section  discusses  the  extension  of  the  SP  model  to  higher 
charge— discharge  rates.  Several  improvements  are  made  to 
different  quantities  including  electrolyte  concentration  and 
potential,  solid  phase  concentration.  Electrodes’  potentials  and  the 
pore  wall  flux  for  higher  rates  are  also  addressed. 

2.3.1.  Electrolyte  concentration 

One  of  the  error-prone  assumptions  of  SP  model  is  neglecting 
the  variation  of  solution  phase  concentration  with  time  and  posi¬ 
tion.  Fig.  2  presents  the  distribution  of  electrolyte  concentration 
inside  the  cell  obtained  by  full  order  model  for  different  discharge 
rates  at  the  end  of  discharge  where  the  cell  voltage  drops  below 


Fig.  2.  Electrolyte  Li  ion  concentration  inside  cell  at  end  of  discharge  for  different  rates, 
P2D  model. 


S.  Khaleghi  Rahimian  et  al.  /  Journal  of  Power  Sources  224  (2013)  180-194 


185 


3.0  V.  The  small  variation  of  concentration  for  lower  rates  (i.e.  0.5C 
and  1.0C)  validates  the  SP  model  constant  concentration  assump¬ 
tion.  However,  as  the  discharge  rate  increases  above  1.0C,  the 
concentration  in  the  electrolyte  changes  remarkably  as  shown  in 
Fig.  2.  Thus,  one  has  to  include  the  solution  phase  material  balance 
(Eq.  (3))  in  SP  model.  The  electrolyte  concentration  profiles  can  be 
approximated  by  polynomial  functions  [13]: 


bruggp 


9Ce,p  (Xp,  f) 


0Xn 


Xp - 1 


fsrUggs9Ce,s(Xs,t) 

Ls  dxs  Xs=o 


£bruggs  ace  s(xs,  t) 

Ls  0XS  Xs=i 


,rSg"0Ce,n(Xn,f) 


Ln 


dxn 


Xn=0 


(19e) 


(19f) 


Ce,i(x, f)  =  cij(t)xm  +  bj(t)xm  1  4  +Zj(t),  jj  _  12  3  (15) 

Different  polynomial  orders  can  be  selected  for  electrolyte 
concentrations  in  the  positive  electrode,  the  separator  and  the 
negative  electrode.  In  this  work,  a  cubic  polynomial  is  chosen  for 
the  electrolyte  concentrations  inside  the  electrodes  while  the 
separator  electrolyte  concentration  is  approximated  by  a  quadratic 
profile: 


The  last  two  equations  are  a  result  of  the  continuity  of 
concentration  and  temperature  at  the  interfaces  and  consequently 
the  diffusion  coefficient  in  Eqs.  (5d)  and  (5f). 

In  order  to  obtain  the  electrolyte  concentration  profile  inside  the 
cell,  11  equations  are  needed  to  find  the  11  polynomial  coefficients 
of  Eq.  (16).  The  boundary  conditions  (Eq.  (19))  provide  6  algebraic 
equations: 

C\  (t)  =  0  (20a) 


Ce,p(xp,t)  =  a,(t)Xj*  +b1(t)x2  +  ci  (t)Xp  +  di(t)  (16a) 

Ce,s(xs,  t)  =  a2(t)x l  +  b2(t)x  s  +  c2(t)  (16b) 

Ce,n(Xn,f)  =  03(f)X^  +  b3(t)X^  +  C3(t)Xn  +  d3(t)  (16c) 

where  x\  is  dimensionless  length  for  each  region: 
x 

xp  =  —  for  0  <  x  <  Lp 
Lp 

<  xs  =  for  Lp  <  x  <  Lp  +  Ls  (17) 

Ts 

xn  =  — — y3 — —  for  Lp  +  Ls  <  x  <  L 

Ln 

The  material  balances  (Eq.  (3))  and  boundary  conditions  (Eq.  (5)) 
can  be  rewritten  using  the  dimensionless  lengths  (Eq.  (17)): 
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9Ce,p  (Xp ,  t) 
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(18a) 

(18b) 


(18c) 

(19a) 


0Ce,n(Xn,  f)  _  ^ 

0Xn  x„=l 

Ce,p(Xp  =  1 ,  f)  =  Ce,s(Xs  =  0,1) 
Ce,s(Xs  =  1,1)  =  Ce_n(Xn  =  0,  1) 


(19b) 

(19c) 

(19d) 


3a3(t)  +  2b3(t)  +  c3(t)  =  0  (20b) 

ai(t)  +  b1(t)  +  c1(t)  +  d1(t)-c2(t)  =0  (20c) 

02(f)  +  b2(t)  +  c2(t)  -  d3(t)  =  0  (20d) 

eLps(3ai (t)  +  2b,  (t)  +  c,  (t))  -  b2(t)  =  0  (20e) 

£Lsn(2a2(f)  +  b2(t))  -  c3(t)  =  0  (20f) 

where  eLps  and  eLsn  are  given  as: 
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Volume  averaging  of  Eq.  (18)  yields  3  following  ordinary 
differential  equations  (ODEs): 
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Xn — 0 
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£nFSn 


(21) 


where  the  average  concentrations  are  functions  of  the  polynomial 
coefficients: 


Q,p,avg  =  (t)  +  —  b\  (t)  +  —  C\  (t)  +  d\  (t) 

Ce,s,avg  =  ^  a2(t)  +  J ^2(0  +  c2(0  (22) 

Q,n,avg  =  ^  <*3  (0  +  3  ^3  (0  +  2  C3  (0  +  ^3  M 

The  two  remaining  equations  can  be  found  by  evaluation  of  the 
material  balances  inside  the  cathode  (Eq.  (18a))  and  the  anode  (Eq. 
(18c))  at  an  interior  point  along  the  x  direction: 
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+ ^ 

xp  =xj 


dCe;n  _  1  9  aCe,n\ 

dt  Xn=x’  L^endXn  \  Oxn  ) 

where jf  ( t )  (i  =  p,n)  represents  the  pore  wall  flux  at  interior  point  x] . 
The  optimal  location  of  the  interior  point  within  the  electrodes  is 
discussed  in  the  Results  and  discussion  section  of  this  manuscript. 


(23b) 
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xP=xi  Lp£p  axp 
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2.3.2.  Electrolyte  potential 

The  electrolyte  potential  distribution  inside  the  cell  for  different 
discharge  rates  at  the  end  of  discharge  is  depicted  in  Fig.  3  illus¬ 
trating  the  reason  for  SP  model  deviation  from  P2D  model  at  higher 
rates  (>1.0C).  Thus,  consideration  of  the  charge  balance  in  the 
solution  phase  is  essential,  in  order  to  predict  the  cell  voltage  at 
higher  rates.  The  solution  phase  potential  inside  each  electrode  is 
also  approximated  by  polynomials  as  follows: 

0e,p(*p?O  =  (0*p  +  (0*p  +  Q  (0*p  +  (0  (24a) 

0e,n(*n,  0  =  A3(t)x*  +  B3(t)x*  +  C3(t)xn  +  D3(t)  (24b) 


while  the  analytical  solution  of  the  electrolyte  potential  in  sepa¬ 
rator  is  obtained  by  assuming  electrolyte  constant  specific 
conductivity.  The  polynomial  coefficients  of  the  electrolyte  poten¬ 
tial  in  the  anode  (i.e.,  A3,  83,  C3,  D3)  are  determined  using  two 
boundary  conditions  at  anode  current  collector  (Eq.  (25a, b)),  the 
charge  balance  at  x\  (Eq.  (25c))  and  volume  averaging  of  both  side 
of  the  charge  balance  (Eq.  (25d)): 


a0e,n 


5xn 
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9  0e,n 
9x2 


(  1  ace,n\ 
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(25b) 
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Fig.  3.  Electrolyte  potential  distribution  inside  cell  at  end  of  discharge  for  different 
rates,  P2D  model. 
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(25d) 


It  is  also  important  to  note  that  in  ISP  model,  the  electrolyte 
potential  at  the  right  end  of  the  cell  is  set  to  zero  (Eq.  (25b))  while 
the  potential  reference  point  in  P2D  model  is  the  solid  potential  at 
the  right  end  of  the  cell.  Since  the  potential  difference  is  considered 
in  the  BV  kinetic  expression,  changing  the  reference  point  does  not 
affect  the  results.  The  charge  balance  in  separator  (Eq.  (6b))  together 
with  the  boundary  conditions  at  the  anode/separator  interface  (Eqs. 
(8e,f))  can  be  rewritten  as: 


a20e,s  o  a  /  1  9Ce,s\ 
6x|  0xs  \Ce,s  axs  / 


=  0 
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(26c) 

Substituting  the  separator  concentration  polynomial  (Eq.  (16b)) 
in  the  boundary  value  problem  (Eq.  (26)),  </>e s  is  expressed  as: 


0e,s 
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where  bc\  and  bc2  are  the  boundary  conditions  (Eq.  (26b, c)): 
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In  order  to  find  the  cathode  electrolyte  potential  coefficients  (i.e., 
A\,  B\,  Ci  and  Di )  boundary  conditions  at  cathode/separator  interface 
(Eq.  (27a, b)),  boundary  condition  at  the  cathode  current  collector 
(Eq.  (27c))  and  the  charge  balance  atxj  (Eq.  (27d))  are  used: 
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In  deriving  Eq.  (25c)  and  Eq.  (27d)  the  solution  phase  specific 
conductivity  was  assumed  to  be  constant  in  x-domain. 


2.3.3.  Solid  potential 

SP  model  assumes  that  the  electrode  potential,  </>;,  depends  on 
time  only.  This  assumption  is  still  true  for  higher  rates  (e.g.  5C)  as 
shown  in  Fig.  4  for  the  cathode  and  anode.  High  conductivity  of  the 
solid  versus  the  electrolyte  makes  the  solid  potential  variation 
inside  the  cell  is  negligible  with  respect  to  the  potential  change  in 
the  solution  phase. 
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Fig.  4.  Potential  distribution  inside  electrodes  at  5C  discharge,  P2D  model  (a)  cathode 
(b)  anode. 


2.3.4.  Solid  concentration 

Instead  of  using  volume  average  technique  for  approximation  of 
the  Li  ion  concentration  inside  particles,  an  approximate  solution 
developed  by  M.  Guo  and  R.E.  White  [18]  is  applied  in  ISP  model.  The 
analytical  solution  for  the  solid  concentration  has  the  following 
form: 


Xs(r,x,  t) 


=  Xs,avg(X,  T)  +  ^5(X,  T)  - 


sin(Anf) 

n4^rsin(A„) 


Qn(X,  t)  + 


2d(x,  t) 


where  xs,aVg(x.T)  and  Qn(x,i)  are  calculated  by  integrating  the 
following  ODEs: 
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where  the  dimensionless  variables  are: 
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It  was  shown  in  Ref.  [18]  that  only  a  few  terms  are  needed  in  the 
analytical  solution  together  with  an  additional  term  compensating 
for  the  truncation  error: 


N 

Xs.surfto'O  =  *s,avg(*,T)  +  £  Qn(*,  l)  +  eN(X,  z) 
n  =  1 

<  eiv(x, t);  =  — 25(x, t)  ^ _  XJ^(1-exp(-Aw+1T)) 
+Y^erfc(AN+1yt)^ 

(29) 

Surface  concentration  can  be  expressed  explicitly  as  a  function 
of  initial  average  concentration  for  constant  current  charge  and 
discharge  where  the  dimensionless  flux  is  constant  (<5(x,t)  =  5).  In 
this  work,  only  one  term  of  approximate  solution  is  applied  to 
calculate  the  surface  concentration: 


Xs,surf,iM  =  x°avg,i  -  35z  -  25  f  (1  -  exp (~xl%) ) 

+y^erfc(A2^)^  -  p  (l  -  exp(-A^)), 

i  =  P,  n 

(30) 

While  the  total  pore  wall  flux  for  each  electrode  (Jj)  is  constant 
during  the  constant  current  charge-discharge,  the  pore  wall  flux  at 
electrodes’  interior  points  (Jj)  vary  with  time.  Thus,  Eq.  (29)  is  used 
to  obtain  the  surface  concentration  by  substituting  xSiavg  and  Qn 
calculated  by  integrating  Eq.  (28)  over  time.  As  a  result,  if  only  one 
term  of  approximate  solution  is  considered,  two  ODEs  (one  for 
average  concentration  and  one  for  Qi )  for  each  surface  concentration 
are  needed. 


2.3.5.  Pore  wall  flux 

In  the  SP  model,  the  applied  current  to  each  particle  (electrode) 
is  equal  to  the  total  current  density  over  the  total  electrode  active 
surface  area  for  that  electrode.  Thus,  the  variation  of  the  pore  wall 
flux  along  x  axis  for  each  electrode  is  neglected  causing  the  SP 
model  cell  voltage  prediction  to  deviate  considerably  from  the 
rigorous  model  at  higher  rates.  Figs.  5  and  6  present  the  variation  of 
the  pore  wall  flux  and  current  flow  in  the  solution  phase,  respec¬ 
tively,  obtained  by  the  full  order  model  inside  the  electrodes  at  a  5C 
discharge.  In  most  reduced  models  [13,16,17]  the  uniform  pore  wall 
flux  is  considered  for  each  electrode  resulting  in  linear  variation  of 
the  current  flow  inside  the  electrode  between  zero  to  the  applied 
current.  Fig.  6  shows  the  deviation  of  the  current  flow  from 
a  straight  line  especially  at  the  end  of  discharge.  In  order  to  increase 
the  accuracy  of  the  model  prediction,  one  has  to  apply  higher  order 
polynomials  for  approximation  of  electrolyte  profiles.  The  BV 
kinetic  can  be  used  at  some  optimal  interior  points  e.g.  (xJ,Xp)  to 
find  the  polynomials’  coefficients: 
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Fig.  5.  Cathode  pore  wall  flux  (a)  and  current  flow  (b)  at  5C  discharge,  P2D  model. 
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Fig.  6.  Anode  pore  wall  flux  (a)  and  current  flow  (b)  at  5C  discharge,  P2D  model. 
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2.3.6.  DAE  system 

Substituting  the  electrolyte  average  concentrations  (Eq.  (22)), 
polynomial  approximation  (Eq.  (16)),  effective  diffusion  (Eq.  (4)) 
and  Eq.  (20a)  in  material  balances  (Eq.  (21))  gives  the  following 
ODEs  in  terms  of  polynomial  coefficients: 


where  xjsurfi,  CF  and  0^-  denote  solid  surface  concentration, 
solution  concentration  and  electrolyte  potential  at  interior  point  xj . 
In  order  to  find  the  electrodes’  potentials,  similar  to  the  SP  model, 
the  BV  expression  is  also  considered  for  the  total  pore  wall  flux  for 
each  electrode: 


where  xStSurftj  is  calculated  using  Eq.  (30)  and  Ce,j,avg  and  0e,i, avg  are 
the  volume  average  electrolyte  concentration  and  potential  inside 
each  electrode  respectively. 
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Eq.  (23)  is  rewritten  using  the  polynomial  approximation  (Eq. 
(16)),  effective  diffusion  (Eq.  (4)): 


vi3dai(f)  udbt(t) 

Xp  dt  +Xp  dt 


d  dj(t) 
dt 


brugg  -1 

fpFSp 


l3da3(t)  l2db3(t)  idc3(t)  d  d3(t) 


dt 


-  +  x 


dt 


+  x. 


bruggn-l 


L2 


Sn(xJ,t) 


+ 


n  dt  +  dt 
(1  -  t°)an/n(t) 


^nFSn 


where 


(34) 


gp(xp,t)  =  aPpd^p’^  (3a1(t)xg+2b1(t)xp) 


-  Dp  (Xp,  t)  (6di  (t)xp  +  2b-[  (t)) 


'ci(t)  =  0 

b2(t)  =  cLps(3a1(t)  +  2ft1(t)  +  c1(t)) 
c2(t)  =  Oi  (t)  +  t>i(t)  +  Ci(t)  +  di  (t) 

<  c3(t)  =  eLsn(2a2(t)  +  b2(t)) 
d3(t)  =  a2(t)  +  b2(t)  +  c2(t) 
b3(t)  =  -(3a3(0  +  c3(t)) 

Thus,  the  system  of  DAEs  includes  5  ODEs  for  5  electrolyte 
polynomial  coefficients,  4  ODEs  (Eq.  (28))  for  solid  average 
concentration  and  Ch*,  2  algebraic  expressions  (Eq.  (32))  for  pore 
wall  fluxes  at  the  interior  points  (J^  Jj)  and  2  algebraic  equations 
(Eq.  (31))  for  solid  phase  potentials.  Thus,  the  new  reduced-order 
model  consists  of  13  DAEs  which  are  solved  in  COMSOL  Inc.  Mul¬ 
tiphysics  3.5a  using  Global  Equations  option  to  compare  the 
computation  time  with  P2D  model.  In  order  to  enhance  the  accu¬ 
racy  of  ISP  model,  higher  polynomial  order  for  the  solution  phase 
concentrations  can  be  used  where  the  evaluation  of  the  material 
balances  at  additional  interior  points  are  required  to  find  the 
polynomial  coefficients.  Since  the  BV  kinetic  expression  is  evalu¬ 
ated  at  additional  points  (Eq.  (31)),  we  are  also  able  to  increase  the 


gn(Xn,t)  =  ^^^-(3a3(t)xl  +  2b3(t)Xn  +  C3(t)') 

+  Dn(xn,  t)(6a3(f)xn  +  2  b3(t)) 

The  derivative  of  De(Xi,t)  with  respect  to  x2-  is  determined  using 
the  chain  rule: 


aPjfo-.t) 

3x, 


/ 8Cej (Xj ,  t)\  ( _ -2.7  X  105 _ 

V  dXj  J  ((7  _  229  -  5.0  x  103Cei(Xj,  t))2 

-0.22  x  103  )(lnlO)Dj(X;,t) 


Taking  the  time  derivative  of  Eq.  (20)  results  in: 


dCT(t)  =  0 
dt 

^da3(t)  +  2db3(t)  i  dc3(t) 


dt 

da-i(t) 

~dT 

da2(t) 

~dT 


dt 


+  - 


=  0 


dfri(t) 

dt 

db2(t) 

dt 


dt 

dci(t)  dd^t) 


el 


VS 


3d^0  +  2 


dt 

dc2(t) 

d  MO 


dt 


dt 

dd3(t) 

dt 

dci(t) 

■”d r 


dc2(t) 

dt 


=  0 


(35) 


=  0 


db2(t) 

dt 


m  o 


da2(t)  ,  d b2(t)\  dc3(t) 


dt 


dt 


dt 


=  0 


Eqs.  (33)— (35)  are  used  to  find  the  right  hand  side  of  the  poly¬ 
nomial  coefficients’  time  derivatives.  At  least  5  ODEs  are  required  to 
integrate  during  charge-discharge  time  for  5  polynomial  coeffi¬ 
cients  (e.g.  ai(t),  bi(t),  di(t),  a2(t),  a3(t)).  The  remaining  6  coeffi¬ 
cients  are  evaluated  in  terms  of  the  5  coefficients  using  Eq.  (20)  as 
follows: 


Fig.  7.  Cell  voltage  maximum  absolute  error  versus  different  interior  points  inside 
electrodes,  (a)  charge  (b)  discharge. 
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Fig.  8.  Cell  voltage  at  different  charge  rates,  (a)  comparison  of  SP  and  ISP  models  with 
P2D  model  (b)  ISP  model  relative  error  with  respect  to  P2D  model. 


Fig.  9.  Cell  voltage  at  different  discharge  rates,  (a)  comparison  of  SP  and  ISP  models 
with  P2D  model  (b)  ISP  model  relative  error  with  respect  to  P2D  model. 


polynomial  order  of  the  electrolyte  potential  as  well  as  Li  ion 
concentrations  in  solution  phase. 

3.  Results  and  discussion 

This  section  compares  the  results  of  original  SP  and  ISP  models 
with  P2D  model  for  different  charge-discharge  rates.  DAEs  systems 
of  the  models  are  solved  by  COMSOL  Inc.  Multiphysics  3.5a  on 
Windows  operating  system  and  the  PC  is  a  Dell  Precision  T7500, 
with  2  Quad  Core  2.53  GHz  Zenon  Processors  CPUs  and  12.285  GB  of 
RAM.  Note  that  only  one  CPU  core  is  involved  in  numerical  simu¬ 
lations  (COMSOL  processor  affinity  is  set  to  one  CPU  (e.g.  CPU  0)  in 
Windows  Task  Manager).  Optimal  interior  points  introduced  in  Eq. 
(23)  are  found  by  applying  MATLAB®  Global  Optimization  Toolbox 
by  minimizing  the  maximum  absolute  error  between  ISP  and  P2D 
models.  The  optimal  points  during  charge  for  all  rates  (i.e.,  1C-5C) 
are  obtained  as  xj  =  0.68  and  x\  =  0.23  while  for  all  rates  (i.e., 
1 C-5C)  during  discharge  the  best  interior  points  are  xj  =  0.06  and 
x^  =  0.12.  The  difference  in  location  of  the  optimal  points  during 
charge  and  discharge  is  due  to  the  difference  in  the  electrolyte 
concentration  and  potential  profiles  between  charge  and  discharge 
throughout  the  cell.  In  order  to  show  the  model  sensitivity  with 


different  interior  points,  the  cell  voltage  error  for  all  rates  (i.e.,  1C— 
5C)  is  plotted  in  Fig.  7  versus  different  interior  points  during  charge 
and  discharge.  The  model  output  does  not  change  considerably 
when  xj  is  between  0.01  and  0.29  or  0.62-0.99  during  charge 
(Fig.  7a)  (0.017  V  <  voltage  maximum  absolute  error  <  0.021  V)  and 
between  0.02  and  0.39  during  discharge  (Fig.  7b)  (0.02  V  <  voltage 


Table  2 

Model  results  for  different  charge  rates,  (a)  SP  and  ISP  models  voltage  errors  with 
respect  to  P2D  model,  (b)  comparison  of  ISP  and  P2D  model  computation  time. 


(a) 

Model 

Voltage  sum  of 

Voltage  maximum 

squared  errors 

absolute  error  (V) 

SP 

3.428 

0.1995 

ISP 

0.019 

0.0169 

(b) 

Model 

CPU  time  (s) 

0.5C  1C  2C 

3C 

4C 

5C 

Mean 

ISP 

2.7  2.3  2.3 

1.92 

2.1 

2.1 

2.2 

P2D 

11.1  10.8  10.1 

12.6 

13.4 

9.2 

11.2 
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Table  3 

Model  results  for  different  discharge  rates,  (a)  SP  and  ISP  models  voltage  errors  with 
respect  to  P2D  model,  (b)  comparison  of  ISP  and  P2D  model  computation  time. 


(a) 

Model 

Voltage  sum  of 
squared  errors 

Voltage  maximum 
absolute  error  (V) 

SP 

8.912 

0.4316 

ISP 

0.016 

0.0198 

(b) 

Model 

CPU  time  (s) 

0.5C  1C  2C 

3C 

4C 

5C 

Mean 

ISP 

2.1  2.1  1.9 

1.8 

2.0 

2.5 

2.1 

P2D 

10.0  9.85  10.9 

12.5 

15.6 

14.3 

12.2 

maximum  absolute  error  <  0.037  V)  when  x\  is  equal  to  its  optimal 
value.  However,  the  sensitivity  of  the  ISP  model  with  respect  to  x\  is 
remarkable  specifically  during  discharge  as  shown  in  Fig.  7b  (the 
error  is  less  than  0.05  V  only  at  optimal  points,  x\  =  0.12).  Thus, 
the  interior  points  could  be  regarded  as  additional  parameters  in 
the  new  model  where  one  can  apply  narrow  bounds  (e.g.  [0.01  0.3]) 
for  parameter  estimation.  The  thickness  (length)  of  the  cathode, 


Cell  Dimensionless  Length 


Fig.  10.  Comparison  electrolyte  concentration  obtained  by  ISP  and  P2D  models  during 
discharge  (a)  1C  (b)  5C. 


separator,  and  anode  has  been  made  dimensionless  in  the  model 
and  thus  the  position  of  the  optimal  internal  points  is  reported  as 
a  dimensionless  parameter.  However  for  a  different  cell  capacity, 
chemistry,  or  charge/discharge  profile  the  location  of  the  internal 
points  should  be  optimized.  Note  that  in  some  regions  of  Fig.  7  the 
errors  are  not  shown  where  either  the  model  DAEs  do  not  converge 
or  the  error  is  outside  the  plot  range. 

3.1.  Cell  voltage 

Figs.  8  and  9  compare  the  cell  voltage  predicted  by  ISP  model 
with  P2D  and  SP  models  for  different  charge  and  discharge  rates 
respectively.  At  lower  rates  (<1.0C)  SP  model  prediction  agrees  well 
with  the  rigorous  model  while  the  voltage  of  SP  model  deviates 
significantly  from  the  rigorous  model  voltage  as  the  charge- 
discharge  rates  increase  above  1C.  On  the  contrary,  the  ISP  model 
is  in  good  agreement  with  the  full  order  model  at  all  rates. 
Tables  2  and  3  summarize  the  results  of  the  ISP  model  compared 
with  the  SP  and  P2D  models  for  charge  and  discharge,  respectively. 
The  residuals  between  the  new  reduced  model  and  P2D  model  are 
much  less  than  the  original  SP  model  errors.  On  the  other  hand,  the 
reduction  in  computational  burden  of  the  rigourous  model  by  using 
the  proposed  model  is  significant. 


Fig.  11.  Comparison  electrolyte  potential  obtained  by  ISP  and  P2D  models  during 
discharge  (a)  1C  (b)  5C. 
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3.2.  Electrolyte  concentration 

Li  ion  concentration  in  the  solution  phase  inside  the  cell  sand¬ 
wich  predicted  by  ISP  model  for  1C  and  5C  discharge  rates  is 
compared  with  electrolyte  concentration  obtained  by  P2D  model  in 
Fig.  10.  Although,  the  new  model  prediction  of  the  electrolyte 
concentration  at  1C  agrees  well  with  the  rigorous  model,  the 
prediction  of  ISP  model  at  5  C  differs  from  the  P2D  model  result. 
However,  the  errors  in  the  solution  phase  concentration  of  the 
reduced  model  do  not  affect  the  cell  voltage  significantly  as  shown 
in  Fig.  9.  In  order  to  improve  the  accuracy  of  ISP  model  electrolyte 
concentration,  higher  order  polynomials  can  be  applied. 

3.3.  Electrolyte  potential 

The  comparison  of  solution  phase  potential  inside  the  cell  ob¬ 
tained  by  ISP  and  P2D  models  at  1C  and  5C  discharge  rates  is 
depicted  in  Fig.  11.  At  1C,  ISP  model  prediction  agrees  well  with  P2D 
model  results  (<0.001  V).  However,  the  difference  between  the  two 
models’  potentials  rises  at  5C  (~0.05  V)  as  shown  in  Fig.  lib.  In 
spite  of  ISP  model  electrolyte  potential  deviation  from  the  full  order 


Fig.  12.  Comparison  of  cell  voltage  at  5C  obtained  by  P2D  and  ISP  models  where  new 
approximation  and  volume  average  technique  are  used  for  solid  concentration  (a) 
charge  (b)  discharge. 


model,  the  cell  voltage  prediction  of  both  model  are  in  good 
agreement  (<1%). 

3.4.  Solid  concentration 


Fig.  12  compares  the  P2D  cell  voltage  at  5C  charge-discharge 
rate  with  the  reduced  model  where  the  volume  average  tech¬ 
nique  (dashed  line)  and  the  new  approximation  [18]  (solid  line)  are 
used  to  predict  the  solid  surface  concentration.  The  improvement 
of  the  cell  voltage  prediction  by  applying  the  new  approximation 
[18]  instead  of  two  term  polynomial  approximation  is  obvious 
especially  at  the  beginning  of  charge  (Fig.  12a)  and  the  end  of 
discharge  (Fig.  12b). 


3.5.  Pore  wall  flux 


The  total  flux  of  Li  ionsj?  (i  =  p,n),  does  not  change  during 
constant  current  charge  and  discharge  (e.g.  for  1C, 

[Jl\  =  ^  =  1.5  x  10-5  rT1°'.  \jl\='f^  =  23  x  10-^). 
rpl  FSP  m2s  r n I  FSn  m2s 

Nonetheless,  the  pore  wall  fluxes  obtained  by  P2D  and  ISP  models 


at  interior  points  inside  the  electrodes  (/J  Jj)  vary  with  time  as 
presented  in  Figs.  13  and  14  for  charge  and  discharge,  respectively. 
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Fig.  14.  Comparison  of  pore  wall  flux  absolute  values  at  interior  points  (xj  =  0.06, 
xj  =  0.12)  obtained  by  ISP  and  P2D  models  during  discharge  (a)  1C  (b)  5C. 

The  ISP  model  fluxes  at  1C  rate  agrees  well  with  the  full  order 
model  fluxes  as  shown  in  Figs.  13a  and  14a  for  charge  and 
discharge,  respectively.  However,  the  difference  between  the  fluxes 
of  two  models  is  considerable  at  5C  specifically  for  the  negative 
electrode  during  charge  (Fig.  13b).  The  pore  wall  flux  error  at  the 
anode  interior  point  can  be  described  by  Figs.  10b  and  lib,  where 
the  ISP  model  electrolyte  concentration  and  potential  differs  from 
the  P2D  model  results.  In  order  to  improve  the  ISP  model  predic¬ 
tion  of  the  interior  pore  wall  flux  inside  the  negative  electrode 
higher  order  polynomial  functions  are  required  for  the  electrolyte 
concentration  and  potential  (e.g.  4th  order  polynomial  for  Ce,n  and 
0e,n)-  Note  also  that  while  the  distribution  of  the  electrolyte 
concentrations  and  potentials  inside  the  cell  can  be  provided  by 
the  ISP  model,  the  pore  wall  flux  can  be  obtained  only  at  the 
interior  points.  Another  point  is  that  we  also  approximated  the 
pore  wall  fluxes  inside  the  electrodes  by  using  polynomials 
profiles.  Then,  the  solution  phase  potentials  were  obtained 
analytically  by  substituting  the  fluxes  approximations  in  the 
electrolyte  charge  balances.  However,  the  model  results  for  elec¬ 
trolyte  quantities  deviated  considerably  with  P2D  model.  Thus, 
polynomial  approximation  of  the  electrolyte  potentials  lead  to 
more  accurate  reduced  model  rather  than  approximation  of  the 
pore  wall  fluxes. 


4.  Conclusion 

An  electrochemical  based  reduced  order  model  has  been 
proposed  by  improvement  of  the  SP  model.  Neglecting  solution 
phase  electrochemical  variations  is  the  main  error  prone  assump¬ 
tion  of  the  SP  model,  resulting  in  significant  error  (>10%)  in  cell 
voltage  prediction  at  higher  rates  (>1C).  In  the  ISP  model  presented 
here,  polynomial  profiles  are  applied  to  approximate  the  Li  ion 
concentrations  and  potential  in  the  solution  phase.  The  BV  kinetic 
expression  in  the  ISP  model  is  both  used  in  average  form  (like  the 
original  SP  model)  as  well  as  at  the  optimal  interior  points  for  each 
electrode.  The  average  concentrations  and  potentials  are  used  in  the 
former  while  the  latter  applies  the  concentrations  and  potentials  at 
the  interior  points.  Hence,  pore  wall  fluxes  inside  the  electrodes  are 
no  longer  uniform  in  the  new  reduced  ISP  model.  Using  more 
accurate  approximation  for  the  solid  phase  diffusion  rather  than 
applying  volume  average  technique,  is  another  enhancement  in  ISP 
model.  Comparison  results  shows  that  ISP  model  is  able  to  predict 
the  cell  voltage  less  than  1%  for  different  charge-discharge  rates  up 
to  5C,  while  the  average  computation  cost  is  saved  by  a  factor  of  5. 
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List  of  symbols 

ai:  specific  interfacial  area  [1/m] 

bruggi:  Bruggeman  number 

Cey.  electrolyte  concentration  [mol  m-3] 

Csf  solid  phase  concentration  [mol  m-3] 

Q.max.i-  maximum  solid  phase  concentration  [mol  m-3] 

Dsf  solid  phase  diffusion  coefficient  of  Li+  [m2  s-1] 

Def.  electrolyte  phase  diffusion  coefficient  of  Li+  [m2  s-1] 

F:  Faraday’s  constant  [C  mol-1] 

Iapp:  applied  current  [C  s-1] 

J,-;  pore  wall  flux  [mol  m-2  s-1] 

kef.  electrolyte  specific  conductivity  [S  m-1] 

lq:  reaction  rate  constant  [m2  5  mol-05  s-1] 

I:  thickness  of  battery  component  [m] 
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polynomial  order 
:  variable  defined  in  Eq.  (28b) 
i .  spatial  coordinate 
r:  dimensionless  radius 
Rg:  gas  constant  [J  mol_1K_1] 

Rf.  particle  radius  (i  =  p,  n)  [m] 

Sf.  electrode  total  active  surface  area  [m2] 
t:  time  [s] 

t°:  transference  number  of  Li  ions  in  electrolyte 
T:  temperature  [K] 

(J,:  open  circuit  potential  [V] 

Vceii:  cell  voltage  [V] 

Vf.  electrode  volume  [m3] 
x:  spatial  coordinate 

*s,avg,i'-  dimensionless  solid  average  concentration 
Xs,surf,i'-  dimensionless  solid  surface  concentration 
fi:  parameter  defined  in  Eq.  (6)  (0  =  2RST(\  -  )/F) 

< Teff,i :  effective  solid  conductivity  [S  m-1] 
rji:  electrode  overpotential  [V] 

</>,■;  electrode  potential  [V] 


{ frej :  electrolyte  potential  [V] 
ei :  volume  fraction  of  active  material 
efi.  volume  fraction  of  fillers 
eLps:  parameter  defined  in  Eq.  (20) 
eLsn:  parameter  defined  in  Eq.  (20) 
An:  eigenvalue 
5:  dimensionless  flux 
t:  dimensionless  time 

Subscripts 
0:  initial  condition 
eff:  efficient 
max:  maximum 
n:  negative  electrode 
p;  positive  electrode 
s:  separator 

Superscripts 

1:  first  interior  point 

T:  total 


